Image processing device, method and program

ABSTRACT

Candidate points belonging to a predetermined structure are extracted from image data DV. A shape model which represents a known shape of the predetermined structure and is formed by model labels having a predetermined connection relationship is obtained. Corresponding points corresponding to the model labels are selected from the extracted candidate points under the following constraints: (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to an image processing device, an image processing method and an image processing program for extracting a plurality of candidate points belonging to a predetermined structure from image data and performing matching between the candidate points and nodes of a known shape model which represents a shape of the predetermined structure.

2. Description of the Related Art

Automatic extraction of a structure, such as coronary arteries, from three-dimensional image data, such as volume data, has conventionally been conducted (see, for example, U.S. Patent Application Publication No. 20110085701, which will hereinafter be referred to as Patent Document 1). Patent Document 1 discloses extracting a plurality of candidate points forming coronary arteries based on values of voxel data forming volume data, selecting, from the extracted candidate points, a corresponding point corresponding to each model label forming a shape model such that the corresponding points form a graph structure that is most similar to the shape model, and determining paths connecting the selected corresponding points by selectively connecting the candidate points to achieve a minimal cost based on a predetermined index value.

The selection of the corresponding points is achieved by providing an evaluation function for evaluating a similarity between the graph structure formed by the set of candidate points corresponding to the model labels and the shape model, and finding an optimal solution (a mapping between the candidate points and the model labels) of the evaluation function. In this case, a set C of feasible solutions x of the evaluation function is usually defined as Equation (1) below:

$\begin{matrix} {C = {\begin{Bmatrix} {{{x \in \left\{ {0,1} \right\}^{P \times Q}}{{\sum\limits_{{Sp} \in P}^{\;}\; {xSpTq}} \leq 1}},{{\sum\limits_{{Tq} \in Q}^{\;}\; {xSpTq}} \leq 1},{\forall{{Sp} \in P}}} \\ {{and}\mspace{14mu} {\forall{{Tq} \in Q}}} \end{Bmatrix}.}} & {(1),} \end{matrix}$

where P represents a set of candidate points S_(p), Q represents a set of model labels T_(q), x_(SpTq) is a binary variable indicating a mapping relationship between an arbitrary candidate point S_(p) and an arbitrary model label T_(q) and has a value of 1 if the candidate point S_(p) is mapped with the model label T_(q) or a value of 0 if the candidate point S_(p) is not mapped with the model label T_(q), and x is a P×Q-dimensional vector with values of the variable x_(SpTq) being the elements thereof.

In Equation (1), the constraints limiting the feasible solutions x are formed by two inequality expressions. The inequality expression on the left means that a sum of values of the variable x_(SpTq) obtained between an arbitrary model label T_(q) and each candidate point S_(p) belonging to the set P is not more than 1. This means that each model label is mapped with only one of the candidate points or none of the candidate points. The inequality expression on the right means that a sum of values of the variable x_(SpTq) obtained between an arbitrary candidate point S_(p) and each model label T_(q) belonging to the set Q is not more than 1. This means that each candidate point is mapped with only one of the model labels or none of the model labels.

However, with the above-described conventional method, where the set of candidate points to be mapped with the model labels is selected from the feasible solutions x belonging to the set C, the selection of the candidate points to be mapped with the model labels as the optimal solution of the evaluation function may be inappropriate, and this is problematic.

For example, as shown in FIG. 11A, when a plurality of candidate points S₁₁ to S₁₉ have been extracted from image data, the candidate points S₁₁, S₁₅, S₁₇ and S₁₉ may be selected as the corresponding points forming a graph structure that matches or is most similar to a shape model M_(ref) 1 formed by model labels T₁₁ to T₁₄, and a graph structure M_(s) 1 may be formed by connecting the corresponding points. However, in the graph structure M_(s) 1, the candidate point selected to be mapped with the model label T₁₂, which is located at the branch point of the structure, is inappropriate, resulting in partially overlapped paths connecting the corresponding points. In this case, it is desired to form a graph structure M_(s) 1′, as shown in FIG. 11B, as an appropriate mapping.

Further, as shown in FIG. 12A, for example, when a plurality of candidate points S₂₁ to S₂₇ have been extracted from image data, the candidate points S₂₁, S₂₆ and S₂₇ may be selected as the corresponding points forming a graph structure that matches or is most similar to a shape model M_(ref) 2 formed by model labels T₂₁ to T₂₃, and a graph structure M_(s) 2 may be formed by connecting the corresponding points. However, in the graph structure M_(s) 2, the candidate point selected to be mapped with the model label T₂₁, which is located at a branch point of the structure, is inappropriate, resulting in overlapped paths connecting the corresponding points. In this case, it is desired to form a graph structure M_(s) 2′, as shown in FIG. 12B, as an appropriate mapping.

SUMMARY OF THE INVENTION

In view of the above-described circumstances, the present invention is directed to providing an image processing device, an image processing method and an image processing program that can achieve more accurate mapping between candidate points extracted from image data and nodes of a shape model.

An aspect of the image processing device of the invention includes: candidate point extracting means for extracting a plurality of candidate points belonging to a predetermined structure from image data; shape model storing means storing a shape model representing a known shape of the predetermined structure, the shape model being formed by a plurality of model labels having a predetermined connection relationship; and corresponding point selecting means for selecting corresponding points corresponding to the model labels from the candidate points extracted by the candidate point extracting means under constraints (a) to (c) below: (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths.

An aspect of the image processing method of the invention includes the steps of: extracting a plurality of candidate points belonging to a predetermined structure from image data; obtaining a shape model stored in advance in shape model storing means, the shape model representing a known shape of the predetermined structure and being formed by a plurality of model labels having a predetermined connection relationship; and selecting corresponding points corresponding to the model labels from the extracted candidate points under constraints (a) to (c) below: (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths.

An aspect of the image processing program of the invention is a program for causing at least one computer to carry out the above-described image processing method. This program is recorded on a recording medium, such as a CD-ROM or DVD, or in a storage of a server computer or a network storage in a downloadable state to be delivered to the user.

The predetermined structure herein may be any object that can be represented by a graph structure formed by points and lines connecting the points, and may include, in particular, a tubular structure of a human body, such as the trachea, the intestines, the coronary arteries, the cerebral vessels, the pulmonary vessels, the hepatic vessels or the bronchial tubes, or a person's face.

The image data herein may be medical image data taken, for example, with a CT device, MR device, ultrasonic device, PET-CT device, SPECT device, 4D-CT device, OCT device or X-ray imaging device (CR, DR), three-dimensional image data, such as volume data, or image data taken with a digital camera, etc.

The shape model herein represents a known shape of a predetermined structure, or, in other words, a shape to be detected of the predetermined structure. For example, the shape model represents a typical anatomical shape of the predetermined structure, or an anatomical shape of the predetermined structure in a normal state.

The paths in the constraint (c) can be determined by selectively connecting the candidate points such that a minimal cost based on a predetermined index value is achieved.

In the image processing device, method and program of the invention, when a set of the candidate points is represented by P, a set of the model labels is represented by Q, and a mapping relationship between the candidate points and the model labels is represented by a vector x having P×Q elements, the corresponding points may be selected by selecting the mapping relationship from a set C of candidate solutions satisfying a constraint expressed by Equation (2) below:

$\begin{matrix} {C = {\begin{Bmatrix} {{{x \in \left\{ {0,1} \right\}^{P \times Q}}{{\sum\limits_{{Sp} \in P}^{\;}\; {xSpTq}} \leq 1}},{{\sum\limits_{{Tq} \in Q}^{\;}\; {xSpTq}} +}} \\ {{{\sum\limits_{{{({{Sp}^{\prime},{Sp}^{''}})} \in {P \times P}}{{Sp}^{\prime} < {Sp}^{''}}}^{\;}\begin{pmatrix} {{ySpSp}^{\prime}{Sp}^{''} \times {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{\prime}{Tq} \times}}} \\ {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{''}{Tq}}} \end{pmatrix}} \leq 1},} \\ {\forall{p \in {P\mspace{14mu} {and}\mspace{14mu} {\forall{q \in Q}}}}} \end{Bmatrix}.}} & {(2),} \end{matrix}$

where x_(SpTq) (x_(Sp′Tq), x_(Sp″Tq)) is a binary variable indicating a mapping relationship between an arbitrary candidate point S_(p) (S_(p′), S_(p″)) and an arbitrary model label T_(q) and has a value of 1 if the candidate point is mapped with the model label or a value of 0 if the candidate point is not mapped with the model label, and y_(SpSp′Sp″) is a binary variable having a value of 1 if the candidate point S_(p) is included in a path between other two candidate points S_(p′) and S_(p″) belonging to the set P or a value of 0 if the candidate point S_(p) is not included in the path.

In Equation (2), the constraints limiting the feasible solutions x are formed by two inequality expressions. The left side of the inequality expression on the left means a sum of values of the variable x_(SpTq) obtained between an arbitrary model label T_(q) and each candidate point S_(p) belonging to the set P, and the first term of the left side of the inequality expression on the right means a sum of values of the variable x_(SpTq) obtained between an arbitrary candidate point S_(p) and each model label T_(q) belonging to the set Q. Further, Σx_(Sp′Tq), which is the second term of the left side of the inequality expression on the right, means a sum of values of the variable x_(Sp′Tq) obtained between the candidate point S_(p′) and each model label T_(q) belonging to the set Q, and Σx_(Sp″Tq) means a sum of values of the variable x_(Sp″Tq) obtained between the candidate point S_(p″) and each model label T_(q) belonging to the set Q.

The corresponding points may be selected by minimizing an energy function E expressed by Equation (3) below including the vector x as a variable using a Dual Decomposition method:

$\begin{matrix} {{{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}},} & (3) \end{matrix}$

where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab) represents a cost of (a,b)εR×R, which is a combination of the mappings a and b.

The extraction of the candidate points may be achieved by detecting a region having an image feature and/or a structural feature of the predetermined structure from the image data, and extracting the candidate points from the detected region.

The selection of the corresponding points may be achieved by any method as long as the corresponding points that match or are most similar to the shape model are selected from the candidate points. For example, a set of corresponding points forming a graph structure that matches or is most similar to the shape model may be selected using graph matching.

The shape model storing means may store an evaluation function, which is obtained through a learning process using training image data known to represent the predetermined structure and indicates a likelihood of the candidate points representing the shape model. In this case, the selection of the corresponding point is achieved using the evaluation function.

The image processing device of the invention may further include structure detecting means for detecting the predetermined structure from the image data based on the corresponding points selected by the corresponding point selecting means and the determined paths between the corresponding points.

The corresponding point selecting means may determine the paths by selectively connecting the candidate points to achieve a minimal cost based on a predetermined index value.

According to the image processing device, method and program of the invention, a plurality of candidate points belonging to a predetermined structure are extracted from image data, a shape model, which is stored in advance in the shape model storing means, represents a known shape of the predetermined structure and is formed by a plurality of model labels having a predetermined connection relationship, is obtained, and corresponding points corresponding to the model labels are selected from the extracted candidate points under the above-described constraints (a) to (c). Therefore, a mapping between the candidate points and the model labels which has overlapped paths connecting the corresponding points, such as the graph structure M_(s) 1 shown in FIG. 11A or the graph structure M_(s) 2 shown in FIG. 12A, can be excluded, thereby achieving more accurate mapping between the candidate points extracted from the image data and the model labels (nodes) of the shape model than that achieved by the conventional method.

Further, in the case where the mapping is selected from a set C of candidate solutions satisfying a constraint expressed by Equation (2) described above, a mapping between the candidate points and the model labels which has overlapped paths connecting the corresponding points can be excluded, thereby achieving more accurate mapping between the candidate points extracted from the image data and the model labels of the shape model.

Still further, in the case where the corresponding points are selected by minimizing the energy function E expressed by Equation (3) described above, which includes the vector x as a variable, using a Dual Decomposition method, fast and accurate selection of the corresponding points can be achieved.

Yet further, in the case where the predetermined structure is detected from the image data based on the corresponding points selected as described above and the determined paths between the corresponding points, accurate detection of the predetermined structure from the image data can be achieved.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram illustrating the schematic configuration of a structure detection device,

FIG. 2 is a schematic diagram illustrating one example of coronary arteries to be detected as a predetermined structure by a candidate region detecting means,

FIG. 3 is a schematic diagram illustrating how candidate points are extracted from a candidate region by a candidate point extracting means,

FIG. 4 is a schematic diagram illustrating how predetermined corresponding points are selected by a corresponding point selecting means,

FIG. 5A is a schematic diagram illustrating one example of training data used for learning a set model,

FIG. 5B is a schematic diagram illustrating one example of training data used for learning a set model,

FIG. 6 is a schematic diagram illustrating one example of a graph (edges) formed by connecting corresponding points (model labels),

FIG. 7 is a schematic diagram illustrating how the corresponding points are select by a corresponding point selecting means,

FIG. 8 is a flow chart illustrating the flow of a process to detect a structure by the structure detection device,

FIG. 9A is a schematic diagram illustrating one example of training data used for learning a set model,

FIG. 9B is a schematic diagram illustrating one example of training data used for learning a set model,

FIG. 10 is a schematic diagram illustrating one example of a graph (edges) formed by connecting the corresponding points (the model labels),

FIG. 11A is a schematic diagram illustrating one example of a graph structure formed by a conventional method,

FIG. 11B is a schematic diagram illustrating a properly formed graph structure for the graph structure shown in FIG. 11A,

FIG. 12A is a schematic diagram illustrating one example of a graph structure formed by a conventional method, and

FIG. 12B is a schematic diagram illustrating a properly formed graph structure for the graph structure shown in FIG. 12A.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

Hereinafter, one embodiment of the present invention will be described with reference to the drawings. FIG. 1 is a block diagram 1 illustrating the schematic configuration of a structure detection device 1. The configuration of the structure detection device 1 as shown in FIG. 1 is implemented by executing, on a computer, a structure detection program, which is read in an auxiliary storage device of the computer. The structure detection program is distributed with being stored in a storage medium, such as a CD-ROM, or via a network, such as the Internet, to be installed on a computer. The structure detection device 1 shown in FIG. 1 detects a structure M, such as coronary arteries, from image data, and includes a candidate region detecting means 10, a candidate point extracting means 20, a normalizing means 30, a corresponding point selecting means 40 and a structure detecting means 50.

The candidate region detecting means 10 detects a candidate region R_(c) by determining whether or not each pixel in the image data DV forms a part of the predetermined structure M. The image data DV is stored in a data storing means VDB. The image data DV may be formed by a two-dimensional image or three-dimensional volume data generated from a plurality of two-dimensional images taken with an imaging device or a radiation detection device, for example.

The candidate region detecting means 10 detects the candidate region using, for example, a technique disclosed in Japanese Patent Application No. 2009-048679 or a technique disclosed in Japanese Patent Application No. 2009-069895, or any other known technique. As one example, a case where the predetermined structure M is the coronary arteries of the heart, as shown in FIG. 2, and the candidate region R_(c) is detected from volume data is described.

First, the candidate region detecting means 10 calculates positions and major axis directions of candidate points forming central lines of the coronary arteries based on values of the voxel data forming the volume data DV. Alternatively, the candidate region detecting means 10 may calculate a Hessian matrix for the volume data DV, and analyze eigenvalues of the calculated Hessian matrix to calculate positional information and the major axis directions of the candidate points forming the central lines of the coronary arteries. Then, the candidate region detecting means 10 calculates, for voxel data around the candidate points, feature quantities indicating likelihood of the voxel data representing the coronary arteries, and determines whether or not the voxel data represents the coronary artery region based on the calculated feature quantities. The determination based on the feature quantities may be achieved, for example, based on an evaluation function, which is obtained in advance through a machine learning process. With this, image data DV, which is determined to represent the coronary artery region, is extracted as the candidate region R_(c) from the image data.

The candidate point extracting means 20 shown in FIG. 1 extracts a plurality of candidate points S_(p) (p=1 to n, where n is the number of extracted candidate points) from the candidate region R_(c) detected by the candidate region detecting means 10. FIG. 3 is a schematic diagram illustrating how the candidate points S_(p) are extracted from the candidate region R_(c). As shown in FIG. 3, the candidate point extracting means 20 extracts a plurality of (primary) candidate points S_(p0) from the candidate region R_(c), and connects the extracted (primary) candidate points S_(p0). For connecting the candidate points, a minimal spanning tree algorithm may be used to connect all the (primary) candidate points S_(p0) at a minimal cost. The cost may be distances between the (primary) candidate points S_(p0), or the like. In a case where the distances between the (primary) candidate points S_(p0) are used as the cost, such a definition may be made that, the larger the distances, the larger the cost. Then, the candidate point extracting means 20 divides the (primary) candidate points S_(p0) per a predetermined length unit (segmentation), and selects the dividing points, where the candidate points are divided, as (secondary) candidate points S_(p). Then, the candidate point extracting means 20 assigns a label S_(p) to each (secondary) candidate point. It should be noted that, although the case where the candidate point extracting means 20 selects the (secondary) candidate points S_(p) from the (primary) candidate points S_(p0) is described as an example, a graph structure M_(s), which will be described later, may be detected using the (primary) candidate points S_(p0). In this case, the candidate region detecting means 10 functions as the candidate point extracting means 20.

The normalizing means 30 shown in FIG. 1 normalizes coordinates of the candidate points S_(p) extracted by the candidate point extracting means 20 based on predetermined reference positions S_(ref). For example, if the predetermined structure is the coronary arteries, the normalizing means 30 detects the aortic valve, the mitral valve and the apex of the heart, as shown in FIG. 2, to set them as the reference positions S_(ref). Then, the normalizing means 30 expresses the position of each candidate point S_(p) in a new coordinate system, where the apex of the heart is the origin, an axis from the apex of the heart toward the aortic valve is a z-axis, an outer product of a vector from the apex of the heart toward the mitral valve and the z-axis is an x-axis, an outer product of the x-axis and, the z-axis is a y-axis, and a length from the apex of the heart to the aortic valve is 1. Since the aortic valve, the mitral valve and the apex of the heart have small variation due to individual differences, using them as the reference positions S_(ref) for the normalization can minimize variation of the candidate points S_(p) due to individual differences.

The corresponding point selecting means 40 shown in FIG. 1 forms the graph structure M_(s) by selecting, from the candidate points S_(p) normalized by normalizing means 30, corresponding points forming the graph structure M_(s) that matches or is most similar to a shape model M_(ref) stored in a shape model storing means DB, and connecting the selected corresponding points such that the graph structure M_(s) substantially matches the shape model M_(ref). As shown in FIG. 4, for example, it is assumed that a plurality of candidate points S_(p) (p=1 to 8) have been extracted by the candidate point extracting means 20, and the shape model M_(ref) formed by four model labels T_(q) (q=1 to 4) having a predetermined connection relationship has been stored in the shape model storing means DB. In this case, the corresponding point selecting means 40 selects the corresponding points forming a graph structure that matches or is most similar to the shape model M_(ref) (i.e., the point S₂ for the label T₁, the point S₄ for the label T₂, the point S₆ for the label T₃, and the point S₇ for the label T₄) to form the graph structure M_(s).

The corresponding point selecting means 40 selects the corresponding points by using graph matching, and the shape model storing means DB stores an energy function E for graph matching as the shape model M_(ref). Then, the corresponding point selecting means 40 finds an optimal solution which achieves minimization of the energy function E from a set C′ of candidate solutions that satisfy the constraint expressed by Equation (4) below:

$\begin{matrix} {C^{\prime} = {\begin{Bmatrix} {{{x \in \left\{ {0,1} \right\}^{P \times Q}}{{\sum\limits_{{Sp} \in P}^{\;}\; {xSpTq}} \leq 1}},{{\sum\limits_{{Tq} \in Q}^{\;}\; {xSpTq}} +}} \\ {{{\sum\limits_{{{({{Sp}^{\prime},{Sp}^{''}})} \in {P \times P}}{p^{\prime} < p^{''}}}^{\;}\begin{pmatrix} {{ySpSp}^{\prime}{Sp}^{''} \times {\sum\limits_{q \in Q}^{\;}\; {{xSp}^{\prime}{Tq} \times}}} \\ {\sum\limits_{q \in Q}^{\;}\; {{xSp}^{''}{Tq}}} \end{pmatrix}} \leq 1},} \\ {\forall{{Sp} \in {P\mspace{14mu} {and}\mspace{14mu} {\forall{{Tq} \in Q}}}}} \end{Bmatrix}.}} & {(4),} \end{matrix}$

where P represents a set of candidate points S_(p), Q represents a set of model labels T_(q), x_(SpTq) (x_(Sp′Tq), x_(Sp″Tq)) is a binary variable indicating a mapping relationship between an arbitrary candidate point S_(p) (S_(p′), S_(p″)) and an arbitrary model label T_(q) and has a value of 1 if the candidate point is mapped with the model label or a value of 0 if the candidate point is not mapped with the model label, and x is a P×Q-dimensional vector with values of the variable x_(SpTq) being the elements thereof.

In Equation (4), the constraints limiting the feasible solutions x are formed by two inequality expressions. The inequality expression on the left means that a sum of values of the variable x_(SpTq) obtained between an arbitrary model label T_(q) and each candidate point S_(p) belonging to the set P is not more than 1. This means that each model label is mapped with only one of the candidate points or none of the candidate points.

The first term of the left side of the inequality expression on the right is a sum of values of the variable x_(SpTq) obtained between an arbitrary candidate point S_(p) and each model label T_(q) belonging to the set Q, and has a value of 1 if the candidate point S_(p) of interest is mapped with any of the model labels T_(q), or has a value of 0 if the candidate point S_(p) of interest is not mapped with any of the model labels T_(q).

The second term of the left side of the inequality expression on the right includes y_(SpSp′Sp″), Σx_(Sp′Tq) and Σx_(Sp″Tq). The variable y_(SpSp′Sp″) is a binary variable. When a path between an arbitrary candidate point S_(p′) and an arbitrary candidate point S_(p″) is determined by selectively connecting candidate points which are mapped with none of the model labels to achieve a minimal cost based on a predetermined index value, the variable y_(SpSp′Sp″) has a value of 1 if an arbitrary candidate point S_(p) is selected to determine the path between the candidate points S_(p′) and S_(p″), or a value of 0 if the candidate point S_(p) of interest is not selected. The variable Σx_(Sp′Tq) is a sum of values of the variable x_(Sp′Tq) obtained between the candidate point S_(p′) and each model label T_(q) belonging to the set Q, and has a value of 1 if the candidate point S_(p′) is mapped with any of the model labels T_(q), or a value of 0 if the candidate point S_(p′) is not mapped with any of the model labels T_(q). The variable Σx_(Sp″Tq) is a sum of values of the variable x_(Sp″Tq) obtained between the candidate point S_(p″) and each model label T_(q) belonging to the set Q, and has a value of 1 if the candidate point S_(p″) is mapped with any of the model labels T_(q), or a value of 0 if the candidate point S_(p″) is not mapped with any of the model labels T_(q).

The inequality expression on the right requires that the total sum of the first term and the second term is not more than 1, that is, one of the first term and the second term is 1 and the other is 0, or both the first term and the second term are 0. That is, when the first term has a value of 1, a solution where the second term has a value of 1 is excluded. This means that each candidate point S_(p) is mapped with only one of the model labels T_(q) or none of the model labels T_(q), and that, when a path between two candidate points S_(p′) and S_(p″) which are mapped with the model labels is determined by selectively connecting the candidate points which are mapped with none the model labels to achieve a minimal cost based on a predetermined index value, each candidate point S_(p) which is mapped with none of the model labels is included in only one of the determined paths or included in none of the determined paths.

The energy function E is expressed by Equation (5) below, which includes the vector x as a variable:

$\begin{matrix} {{{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}},} & (5) \end{matrix}$

where R represents a set of mappings between the candidate points and the model labels in the feasible solutions x, and θ_(a) is a cost of an arbitrary mapping a (between an arbitrary candidate point S_(p′) and an arbitrary model label T_(q′)) belonging to the set R, which is expressed by Equation (6) below:

θ_(a) =−L _(Tq′)(C _(Sp′))  (6).

In Equation (6), L_(Tq′) (C_(Sp′)) represents a probability of the model label T_(q′) being mapped with the candidate points S_(p′), and the distribution thereof is approximated by a Gaussian function, which is expressed by Equation (7) below:

L _(Tq′)(C _(Sp′))=N(C _(Sp′)|μ_(Tq′),σ_(Tq′) ²)  (7),

where C_(sp′) represents coordinates of the candidate point S_(p′), and μ_(Tq′) and σ_(Tq′) ² are average coordinates and a variance, respectively, of the model label T_(q′), which is obtained from a plurality of pieces of training data known to represent the predetermined structure M. The value of L_(Tq′(C) _(Sp′)) is smaller as a distance between the coordinates of the candidate point S_(p′) and the coordinates of the model label T_(q′) is larger, and vice versa.

It should be noted that the model labels T_(q) (q=1 to 26, for example) are extracted in advance from training data, which is known to represent the predetermined structure M, as shown in FIGS. 5A and 5B. In FIGS. 5A and 5B, the model labels T1 to T20 are extracted from the coronary arteries, which are the predetermined structure M, and the model labels T21 to T26 are extracted from the coronary veins. Then, as shown in FIG. 6, the shape model (graph) M_(ref) is formed by connecting predetermined labels T_(q′) among the model labels T1 to T26 by edges (lines). The extraction of the model labels T_(q) is carried out for a plurality of pieces of training data. Then, for each model label T_(q′) of each training data, average coordinates μ_(Tq′) and a variance σ_(Tq′) ² are calculated, and a probability density function L_(Tq′) (C_(Sp′)) expressed by Equation (7) is derived for each model label T_(q′).

On the other hand, θ_(ab) in Equation (5) is a cost of a combination of two mappings a and b belonging to the set R (a is a mapping between the candidate point S_(p′) and the model label T_(q′), and b is a mapping between the candidate point S_(p″) and the model label T_(q″)), and is expressed by Equation (8) below:

$\begin{matrix} {\theta_{ab} = {{- {L_{{Tq}^{\prime}{Tq}^{''}}\left( C_{{Sp}^{\prime},{Sp}^{''}} \right)}} + {\min {\sum\limits_{{({{Si},{Sj}})} \in {{Path}\mspace{14mu} {({{Sp}^{\prime},{Sp}^{''}})}}}^{\;}{{e_{local}\left( {{Si},{Sj}} \right)}.}}}}} & (8) \end{matrix}$

It should be noted that, in this embodiment, a sum of costs obtained by finding a path connecting each pair of S_(p′) and S_(p″) to achieve a minimal cost is used as the matching cost θ_(ab). However, a different evaluation value relating to the found paths may be used as the matching cost. For example, an evaluation value, such as an average or variance of luminance values of pixels located along the paths, an average or variance of values of the radius of the structure M at each position along the paths, or the entire length of the paths, may be used.

The variance L_(Tq′Tq″) (C_(Sp′)C_(Sp″)) represents a probability of an arbitrary pair of model labels (T_(q′), T_(q″)) being mapped with an arbitrary pair of candidate points (S_(p′), S_(p″)), and is expressed by Equation (9) below:

L _(Tq′Tq″)(C _(Sp′) ,C _(Sp″))=N(C _(Sp′) −C _(Sp″)|μ_(Tq′Tq″),σ_(Tq′Tq″) ²)  (9),

where (C_(Sp′)−C_(Sp″)) is relative coordinates of the candidate point S_(p″) relative to the coordinates of the candidate point S_(p′), and μ_(Tq′Tq″) and σ_(Tq′Tq″) ² are an average value and a variance, respectively, of relative coordinates of the model label T_(q″) relative to the coordinates of the model label T_(q′) obtained from a plurality of pieces of training data which are known to represent the predetermined structure M (see FIG. 7).

Further, e_(local) in the second term of Equation (8) represents a cost of a connection between directly-connected two candidate points S_(i), S_(j), if any, which are located along a path connecting the candidate points S_(p′) and S_(p″) mapped with the model labels, may be defined such that e_(local) has a smaller value as the distance between the two points S_(i), S_(j) is smaller, or the difference of the luminance value between the two points S_(i), S_(j) is smaller, or the difference of the width of the blood vessel between the two points S_(i), S_(j) is smaller, for example. Σe_(local) is a sum of values of the e_(local) obtained for all the pairs of directly-connected two candidate points along the path between the candidate points S_(p′) and S_(p″), and minΣe_(local) is found by calculating an optimal path between the candidate points S_(p′) and S_(p″) using, for example, an Dijkstra's algorithm with handling the candidate points as a graph structure where each side (edge) connecting the candidate points is weighted with a weight defined by the e_(local).

Further, the corresponding point selecting means 40 stores paths connecting the candidate points between the candidate points S_(p′) and S_(p″) which achieve a minimal value of Σe_(local) in a predetermined memory area of the structure detection device 1. With this, an optimal path that smoothly connects the candidate points S_(p′) and S_(p″) via other candidate points is determined from paths connecting various candidate points between the candidate points S_(p′) and S_(p″).

The corresponding point selecting means 40 solves (i.e., finds an optimal solution of) the above-described minimization problem of the constrained energy function E using a Dual Decomposition (DD) method. The DD method divides a main problem to be solved into sub-problems that are calculable. The main problem becomes a convex combination if the condition expressed by Equation (10) below is satisfied:

$\begin{matrix} {{{\sum\limits_{\sigma \in I}^{\;}\; {\rho_{\sigma}\theta^{\sigma}}} = \theta},} & (10) \end{matrix}$

where θ is a vector representing the energy function of the main problem, each sub-problem is represented by a vector θ^(σ) (σεI), ρ_(σ) is a weighting factor having a positive value, and I is an index of the sub-problem.

A lower bound function Φ(θ), which is defined based on the sub-problems, is a convex function expressed by Equation (11) below:

Φ(θ)=Σ_(σεI)ρ_(σ) min E(x|θ ^(σ))  (11),

and can therefore be maximized by a projected subgradient method. There is a relationship that a lower bound value is not larger than a solution (upper bound value) of the main problem, and therefore the minimal solution of the main problem is obtained when the lower bound value is the same value as the upper bound value.

The division into the sub-problems is achieved as follows. For each node, one sub-problem that takes only the node and an adjacent node connected to the node by an edge into account is generated. Each sub-problem ignores the costs θ_(a) and θ_(ab) of the other nodes (i.e., considers the costs as “0”). When the sub-problem is small enough, the solution thereof can be calculated by searching through all the combinations. If a branch-and-bound method is used, the solution can be obtained more rapidly. In the case where the branch-and-bound method is used, first, a relaxed lower bound is set. In this embodiment, the relaxed lower bound is defined as a minimal cost which may possibly be achieved by a certain node without a uniqueness constraint. By omitting detailed calculation of combinations which are not likely to improve the solution with referencing the relaxed lower bound, the amount of calculation can significantly be reduced.

Alternatively, the corresponding point selecting means 40 may solve the minimization problem of the energy function E by using a loopy belief propagation method. Further, although the energy function E is minimized in Equation (5), the first variable term and the second variable term may be multiplied with a negative value, and the evaluation value of the energy function E may be maximized.

The structure detecting means 50 shown in FIG. 1 detects the structure M based on the corresponding points selected by the corresponding point selecting means 40 and the optimal paths between the corresponding points determined in advance. It should be noted that the optimal paths have already been calculated and stored in a predetermined memory area of the structure detection device 1 when the corresponding point selecting means 40 calculates the minΣe_(local).

In a case where the shape model M_(ref) does not represent the entire structure M to be detected, the structure detecting means 50 detects the entire structure M by using a minimal spanning tree algorithm to connect candidate points which are not used in the optimal path between the corresponding points. At this time, when a distance between the candidate points are larger than a predetermined value, the candidate points are not connected, and a group of candidate points which is eventually not connected to any part of the detected structure M detected based on the shape model M_(ref) is handled as not belonging to the structure M. In contrast, in a case where the shape model M_(ref) represents the entire structure M to be detected, it is not necessary to the connect candidate points which are not used in the optimal path between the corresponding points, since they are remaining regions (such as pulmonary vessels or the heart muscle) that do not belong to the structure M to be detected.

FIG. 8 is a flow chart illustrating the flow of a process to detect a structure by the structure detection device. As shown in FIG. 8, first, the candidate region detecting means 10 detects the candidate region R_(c) from the image data DV (step ST1, see FIG. 2). Then, the candidate point extracting means 20 extracts the (primary) candidate points S_(p0) from the candidate region R_(c), and extracts the (secondary) candidate points S_(p) by segmentation (step ST2, see FIG. 3). Thereafter, the normalizing means 30 performs normalization, and the corresponding point selecting means 40 selects, from the candidate points S_(p), the set of corresponding points forming the graph structure M_(s) that is most similar to the shape model M_(ref), according to Equation (5) (step ST3, see FIG. 4). Then, the structure detecting means 50 detects the structure M based on the optimal paths connecting the selected corresponding points via other candidate points (step ST4).

As described above, according to the structure detection device, method and program of this embodiment, the candidate points belong to the predetermined structure M are extracted from the image data DV. From the extracted candidate points, the corresponding points forming the graph structure M_(s) that matches or is most similar to the shape model M_(ref), which is formed by the model labels T_(q) having a predetermined connection relationship, represents a known shape of the structure M and is stored in advance in the shape model storing means DB, are selected under the following constraints (a) to (c): (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths. Then, the structure M is detected from the image data DV using the selected corresponding points. This allows accurate detection of the structure M from the image data DV. In particular, the constraint (c) serves to exclude a mapping between the candidate points and the model labels which has overlapped paths connecting the corresponding points, such as the graph structure M_(s) 1 shown in FIG. 11A or the graph structure M_(s) 2 shown in FIG. 12A, thereby improving the accuracy of the detection of the structure M compared to the conventional method disclosed in Patent Document 1, or the like.

Further, in a case where the corresponding point selecting means 40 selects the corresponding points under the constraint where an arbitrary candidate point S_(p) belonging to set P and an arbitrary model label T_(q) belonging to the set Q satisfy Equation (4) above, a mapping between the candidate point and the model label with overlapped paths connecting the corresponding points can be excluded, thereby allowing accurate detection of the structure M from the image data DV.

Although the above-described embodiment is the case where the image processing device, method and program of the invention are applied to detection of the coronary arteries from the image data DV, the image processing device, method and program of the invention are also applicable to extraction of the trachea, the intestines, the cerebral vessels, the pulmonary vessels, the hepatic vessels or the bronchial tubes, for example. In this case, a shape model corresponding to the structure M, such as the cerebral vessels, is stored in the shape model storing means DB, and the corresponding points are selected using a designed shape depending on the type of the structure M.

For example, in a case where the structure M is the cerebral vessels, the model labels T_(q) (q=1 to 25, for example) are extracted from training data that is known to represent the cerebral vessels, as shown in FIGS. 9A and 9B. In FIGS. 9A and 9B, the model labels T₁ to T₁₂ are extracted from the arteries, which are the predetermined structure M, and the model label T₁₃ to T₂₅ are extracted from the veins. Then, as shown in FIG. 10, the shape model (graph) M_(ref) is formed by connecting predetermined labels among the model labels T1 to T25 by edges (lines).

The image processing device, method and program of the invention are also applicable to a case where graph structures which are extracted from a plurality of images taken by imaging the same patient at different phases are mapped with each other. In this case, one of the graph structures extracted from one of phases may be used as the shape model, and nodes forming the graph structure may be used as the model labels.

Although the above-described embodiment includes the normalizing means 30 and the structure detecting means 50, the normalizing means 30 and the structure detecting means 50 are optional.

In a case where a predetermined structure M including different sub-structures (arteries and veins, for example) is used as the shape model M_(ref) the structure detection device 1 may further include a display control means for displaying the different sub-structures in different colors.

In a case where the predetermined structure is blood vessels, the corresponding point selecting means 40 may select the corresponding points by using coordinates of the candidate points and radii or luminance values of the blood vessels at the candidate points. 

What is claimed is:
 1. An image processing device comprising: a candidate point extracting unit configured to extract a plurality of candidate points belonging to a predetermined structure from image data; a shape model storing unit configured to store a shape model representing a known shape of the predetermined structure, the shape model being formed by a plurality of model labels having a predetermined connection relationship; and a corresponding point selecting unit configured to select corresponding points corresponding to the model labels from the candidate points extracted by the candidate point extracting unit under constraints (a) to (c) below: (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths.
 2. The image processing device as claimed in claim 1, wherein, when a set of the candidate points is represented by P, a set of the model labels is represented by Q, and a mapping relationship between the candidate points and the model labels is represented by a vector x having P×Q elements, the corresponding point selecting unit selects the mapping relationship from a set C of candidate solutions satisfying a constraint expressed by Equation (1) below: $\begin{matrix} {C = {\begin{Bmatrix} {{{x \in \left\{ {0,1} \right\}^{P \times Q}}{{\sum\limits_{{Sp} \in P}^{\;}\; {xSpTq}} \leq 1}},{{\sum\limits_{{Tq} \in Q}^{\;}\; {xSpTq}} +}} \\ {{{\sum\limits_{{{({{Sp}^{\prime},{Sp}^{''}})} \in {P \times P}}{{Sp}^{\prime} < {Sp}^{''}}}^{\;}\begin{pmatrix} {{ySpSp}^{\prime}{Sp}^{''} \times {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{\prime}{Tq} \times}}} \\ {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{''}{Tq}}} \end{pmatrix}} \leq 1},} \\ {\forall{p \in {P\mspace{14mu} {and}\mspace{14mu} {\forall{q \in Q}}}}} \end{Bmatrix}.}} & {(1),} \end{matrix}$ where x_(SpTq) is a binary variable having a value of 1 if a certain one of the candidate points S_(p) is mapped with a certain one of the model labels T_(q) or a value of 0 if the certain one of the candidate points S_(p) is not mapped with the certain one of the model labels T_(q), and y_(SpSp′Sp″) is a binary variable having a value of 1 if the certain candidate point S_(p) is included in a path between other two candidate points S_(p′) and S_(p″) belonging to the set P or a value of 0 if the certain candidate point S_(p) is not included in the path.
 3. The image processing device as claimed in claim 1, wherein the corresponding point selecting unit selects the corresponding points by minimizing an energy function E expressed by Equation (2) below including the vector x as a variable using a Dual Decomposition method: $\begin{matrix} {{{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}},} & (2) \end{matrix}$ where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab) represents a cost of (a,b)εR×R, which is a combination of the mappings a and b.
 4. The image processing device as claimed in claim 1, further comprising a structure detecting unit configured to detect the predetermined structure from the image data based on the corresponding points selected by the corresponding point selecting unit and the determined paths between the corresponding points.
 5. The image processing device as claimed in claim 1, wherein the corresponding point selecting unit determines the paths by selectively connecting the candidate points to achieve a minimal cost based on a predetermined index value.
 6. The image processing device as claimed in claim 2, wherein the corresponding point selecting unit selects the corresponding points by minimizing an energy function E expressed by Equation (2) below including the vector x as a variable using a Dual Decomposition method: $\begin{matrix} {{{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}},} & (2) \end{matrix}$ where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab) represents a cost of (a,b)εR×R, which is a combination of the mappings a and b.
 7. The image processing device as claimed in claim 2, further comprising a structure detecting unit configured to detect the predetermined structure from the image data based on the corresponding points selected by the corresponding point selecting unit and the determined paths between the corresponding points.
 8. The image processing device as claimed in claim 2, wherein the corresponding point selecting unit determines the paths by selectively connecting the candidate points to achieve a minimal cost based on a predetermined index value.
 9. An image processing method comprising the steps of: extracting a plurality of candidate points belonging to a predetermined structure from image data; obtaining a shape model stored in advance in a shape model storing unit, the shape model representing a known shape of the predetermined structure and being formed by a plurality of model labels having a predetermined connection relationship; and selecting corresponding points corresponding to the model labels from the extracted candidate points under constraints (a) to (c) below: (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths.
 10. The image processing method as claimed in claim 9, wherein, when a set of the candidate points is represented by P, a set of the model labels is represented by Q, and a mapping relationship between the candidate points and the model labels is represented by a vector x having P×Q elements, the step of selecting corresponding points comprises selecting the mapping relationship from a set C of candidate solutions satisfying a constraint expressed by Equation (1) below: $\begin{matrix} {C = {\begin{Bmatrix} {{{x \in \left\{ {0,1} \right\}^{P \times Q}}{{\sum\limits_{{Sp} \in P}^{\;}\; {xSpTq}} \leq 1}},{{\sum\limits_{{Tq} \in Q}^{\;}\; {xSpTq}} +}} \\ {{{\sum\limits_{{{({{Sp}^{\prime},{Sp}^{''}})} \in {P \times P}}{{Sp}^{\prime} < {Sp}^{''}}}^{\;}\begin{pmatrix} {{ySpSp}^{\prime}{Sp}^{''} \times {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{\prime}{Tq} \times}}} \\ {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{''}{Tq}}} \end{pmatrix}} \leq 1},} \\ {\forall{p \in {P\mspace{14mu} {and}\mspace{14mu} {\forall{q \in Q}}}}} \end{Bmatrix}.}} & {(1),} \end{matrix}$ where x_(SpTq) is a binary variable having a value of 1 if a certain one of the candidate points S_(p) is mapped with a certain one of the model labels T_(q) or a value of 0 if the certain one of the candidate points S_(p) is not mapped with the certain one of the model labels T_(q), and y_(SpSp′Sp″) is a binary variable having a value of 1 if the certain candidate point S_(p) is included in a path between other two candidate points S_(p′) and S_(p″) belonging to the set P or a value of 0 if the certain candidate point S_(p) is not included in the path.
 11. The image processing method as claimed in claim 9, wherein the step of selecting corresponding points comprises selecting the corresponding points by minimizing an energy function E expressed by Equation (2) below including the vector x as a variable using a Dual Decomposition method: $\begin{matrix} {{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}} & (2) \end{matrix}$ where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab) represents a cost of (a,b)εR×R, which is a combination of the mappings a and b.
 12. The image processing method as claimed in claim 10, wherein the step of selecting corresponding points comprises selecting the corresponding points by minimizing an energy function E expressed by Equation (2) below including the vector x as a variable using a Dual Decomposition method: $\begin{matrix} {{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}} & (2) \end{matrix}$ where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab), represents a cost of (a,b)εR×R, which is a combination of the mappings a and b.
 13. A non-transitory computer readable medium containing an image processing program for causing a computer to carry out the steps of: extracting a plurality of candidate points belonging to a predetermined structure from image data; obtaining a shape model stored in advance in a shape model storing unit, the shape model representing a known shape of the predetermined structure and being formed by a plurality of model labels having a predetermined connection relationship; and selecting corresponding points corresponding to the model labels from the extracted candidate points under constraints (a) to (c) below: (a) each model label is mapped with only one of the candidate points or none of the candidate points; (b) each candidate point is mapped with only one of the model labels or none of the model labels; and (c) when a path between two candidate points which are mapped with each pair of the model labels connected with each other is determined, each candidate point which is mapped with none of the model labels is included in only one of the determined paths or none of the determined paths.
 14. The non-transitory computer readable medium as claimed in claim 13, wherein, when a set of the candidate points is represented by P, a set of the model labels is represented by Q, and a mapping relationship between the candidate points and the model labels is represented by a vector x having P×Q elements, the step of selecting corresponding points comprises selecting the mapping relationship from a set C of candidate solutions satisfying a constraint expressed by Equation (1) below: $\begin{matrix} {C = {\begin{Bmatrix} {{{x \in \left\{ {0,1} \right\}^{P \times Q}}{{\sum\limits_{{Sp} \in P}^{\;}\; {xSpTq}} \leq 1}},{{\sum\limits_{{Tq} \in Q}^{\;}\; {xSpTq}} +}} \\ {{{\sum\limits_{{{({{Sp}^{\prime},{Sp}^{''}})} \in {P \times P}}{{Sp}^{\prime} < {Sp}^{''}}}^{\;}\begin{pmatrix} {{ySpSp}^{\prime}{Sp}^{''} \times {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{\prime}{Tq} \times}}} \\ {\sum\limits_{{Tq} \in Q}^{\;}\; {{xSp}^{''}{Tq}}} \end{pmatrix}} \leq 1},} \\ {\forall{p \in {P\mspace{14mu} {and}\mspace{14mu} {\forall{q \in Q}}}}} \end{Bmatrix}.}} & {(1),} \end{matrix}$ where x_(SpTq) is a binary variable having a value of 1 if a certain one of the candidate points S_(p) is mapped with a certain one of the model labels T_(q) or a value of 0 if the certain one of the candidate points S_(p) is not mapped with the certain one of the model labels T_(q), and y_(SpSp′Sp″) is a binary variable having a value of 1 if the certain candidate point S_(p) is included in a path between other two candidate points S_(p′) and S_(p″) belonging to the set P or a value of 0 if the certain candidate point S_(p) is not included in the path.
 15. The non-transitory computer readable medium as claimed in claim 13, wherein the step of selecting corresponding points comprises selecting the corresponding points by minimizing an energy function E expressed by Equation (2) below including the vector x as a variable using a Dual Decomposition method: $\begin{matrix} {{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}} & (2) \end{matrix}$ where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab) represents a cost of (a,b)εR×R, which is a combination of the mappings a and b.
 16. The non-transitory computer readable medium as claimed in claim 14, wherein the step of selecting corresponding points comprises selecting the corresponding points by minimizing an energy function E expressed by Equation (2) below including the vector x as a variable using a Dual Decomposition method: $\begin{matrix} {{E\left( {x\theta} \right)} = {{\sum\limits_{a \in R}^{\;}{\theta_{a}x_{a}}} + {\sum\limits_{{({a,b})} \in {R \times R}}^{\;}{\theta_{ab}x_{a}x_{b}}}}} & (2) \end{matrix}$ where a and b represent individual mappings between the candidate points and the model labels in the mapping relationship represented by the vector x, b represents a different mapping from a mapping represented by a, a and b are elements belonging to a set R and the set R is a subset of P×Q, θ_(a) represents a cost of a mapping aεR, and θ_(ab) represents a cost of (a,b)εR×R, which is a combination of the mappings a and b. 